package edu.gmu.atelier.generic;

import edu.gmu.atelier.generic.RealVectorUtil;
import edu.gmu.atelier.generic.RealVectorFactory;
import edu.gmu.atelier.generic.RealVector;
import edu.gmu.atelier.generic.AbstractMatrix;
import java.util.ArrayList;
import java.util.List;

/**
 * Performs LU decomposition on the passed matrix.
 * @author James H. Pope
 */
public class LUDecomposition
{
    private AbstractMatrix<RealVector> a = null;
    
    private AbstractMatrix<RealVector> u = null;
    
    private AbstractMatrix<RealVector> l = null;
    
    private int[] perms  = null; // technically permutation matrix
    
    public LUDecomposition( AbstractMatrix<RealVector> a )
    {
        this.a = a;
        this.computeLU();
        //this.l = new RealMatrix(computeL(m));
    }
    
    //------------------------------------------------------------------------//
    // Interface methods
    //------------------------------------------------------------------------//
    public AbstractMatrix<RealVector> getU()
    {
        return u;
    }
    
    public AbstractMatrix<RealVector> getL()
    {
        return l;
    }
    
//    /**
//     * Solves the equation Ax=LUx=b.  Since Ux=y, Ly=b, so we first solve for
//     * y from L and specified b.  Then we solve for x in equation Ux=y.
//     * @return x
//     */
//    public RealVector solve( RealVector b )
//    {
//        RealVector y = LUDecomposition.forwardSubstitution(l, b);
//        //System.out.println( "Solution y: " + y.toString() );
//
//        RealVector x = LUDecomposition.backSubstitution(u, y);
//        //System.out.println( "Solution x: " + x.toString() );
//
//        return x;
//    }
    
//    /**
//     * Gets the determinant of the matrix A (actually determined by using
//     * the already factored upper triangular U and simply multiplying the
//     * diagonal).
//     * @return det(A)
//     */ 
//    public double determinant()
//    {
//        double det;
//
//        // We have already checked that matrix A is square so this is safe.
//        int maxSize = this.u.getRows();
//        det = this.u.cols[0].components[0];
//        for(int i = 1; i<maxSize; i++)
//        {
//            det = det * this.u.cols[i].components[i];
//        }
//
//        return det;
//    }
    
    public String toString()
    {
        StringBuffer buf = new StringBuffer();
        buf.append( "Matrix a: " + a.toString() );
        buf.append( "Matrix u: " + u.toString() );
        buf.append( "Matrix l: " + l.toString() );
        return buf.toString();
    }
    
    private void computeLU( )
    {
        // Make permutation "matrix"
        this.perms = new int[a.getRows()];
        for (int i = 0; i < perms.length; i++)
        {
            perms[i] = i;
        }
        
        // Make u - starting out a copy of a
        List<RealVector> u = new ArrayList<RealVector>();
        for (int i = 0; i < a.getCols(); i++)
        {
            u.add( i, a.getVector(i).copy() );
        }
        
        List<RealVector> l = RealVectorFactory.makeIdentityVectors( u.size() );
        
        for (int i = 0; i < u.size()-1; i++)
        {
            // Partial Pivoting (by row), find largest pivot value in column
            int s = i;
            double largestValue = u.get(i).elements[i];
            for (int z = i; z < u.get(0).getLength(); z++)
            {
                double possiblePivot = u.get(i).elements[z];
                if( possiblePivot > largestValue )
                {
                    largestValue = possiblePivot;
                    s   = z;
                }
            }
            
            //System.out.println("LUDecom "+i+": "+ new AbstractMatrix(u) );
            if(i!=s)
            {
                //System.out.println("Swapping rows "+i+ " and "+s);
                RealVectorUtil.swap(u, i, s);
                // Update the permutation "matrix"
                int ival = this.perms[i];
                this.perms[i] = s;
                this.perms[s] = ival;
            }
            //System.out.println("Swapped  "+i+": "+ new AbstractMatrix(u) );
            
            double pivotCoeff = u.get(i).elements[i];
            for (int r = i+1; r < u.size(); r++)
            {
                double lowerCoeff = u.get(i).elements[r];
                //System.out.println("i="+i+" pivotCoeff="+pivotCoeff+" lowerCoeff="+lowerCoeff);
                double multiplier = lowerCoeff / pivotCoeff;
                RealVectorUtil.mulAdd(0, u, i, -multiplier, r);
                
                l.get(i).elements[r] = multiplier;
            }
        }
        this.u = new AbstractMatrix<RealVector>(u);
        this.l = new AbstractMatrix<RealVector>(l);
        
        System.out.println("Perms:");
        for (int i = 0; i < perms.length; i++)
        {
            System.out.println("      "+i+"="+perms[i]);
        }
        System.out.println("\n");
    }
    
//    private void computeLU( )
//    {
//        double[][] u = RealColMatrixFactory.makeData( m );
//        double[][] l = RealColMatrixFactory.identityData(u.length);
//        for (int i = 0; i < u.length-1; i++)
//        {
//            double pivotCoeff = u[i][i];
//            for (int r = i+1; r < u.length; r++)
//            {
//                double lowerCoeff = u[r][i];
//                //System.out.println("i="+i+" pivotCoeff="+pivotCoeff+" lowerCoeff="+lowerCoeff);
//                double multiplier = lowerCoeff / pivotCoeff;
//                DoubleUtil.mulAdd(0, u, i, -multiplier, r);
//                
//                l[r][i] = multiplier;
//            }
//        }
//        this.u = new RealColMatrix(u);
//        this.l = new RealColMatrix(l);
//    }
    
    
}
